Continuum model of actin-mediated bacterial propulsion 
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The mechanism of propulsion of host bacteria under the action of actin gel networks is examined 
by means of a continuum model of the dynamics of the F-actin concentration. The model includes 
the elasticity of the network, its attachment to the host and the polymerization at the interface with 
it. A formula for the cruise velocity is derived wherefrom the respective contributions of elasticity 
and polymerization are made explicit. The velocity tends to elongate the gel in the direction 
opposite to motion, and in turn is proportional to the size of its active portion. It is shown that the 
motion can start only past a finite latency time: the condition for the onset of motion is explicitly 
given. It is numerically found that at steady state the motion has a pulsating behavior, with sudden 
decelerations and subsequent recoveries. 



I. INTRODUCTION 



The mechanism of actin-induced propulsion has been 
the object of increasingly sophisticated experiments in 
recent years [l|, • The phenomenon is of interest as a 
paradigm of more complex situations, and exhibits a va- 
riety of fascinating and intriguing features, which keep 
stimulating new research. In short, an alien host (e.g. a 
bacterium or a polystyrene sphere) can make the cyto- 
plasmic material to react by polymerizing the resident G- 
actin into an F-actin gel network surrounding the host's 
surface; in due time the network evolves into an asym- 
metric shape, typically comet-like. The host eventually 
starts moving, and its cruise can last for relatively long 
times at almost constant velocity. We will focus here on 
two features of the process; the first one is a threshold 
mechanism at the onset of motion: the velocity of the 
host rises from zero to its cruise value quite abruptly, 
while the mass of the gel undergoes a smooth growth. 
The second feature is a correlation between the cruise 
velocity and the length of the comet [l|:longer tails gen- 
erate higher velocities. This suggests that the dominant 
role in propulsion stems from a bulk force, rather than 
from an interface effect, such as the polymerization rate 
at the host's surface. In fact in cases of growth over con- 
vex surfaces it has been observed and proved that the 
gel elasticity has a primary role: the stress can inhibit 
the polymerization at the interface and become the 
main propulsive factor. Disk-shaped hosts can as well 
be put into motion [4]; one could argue that since the 
gel coating entirely surrounds the host, the convex por- 
tions of its surface are in all cases under the action of 
stress. It is difficult though to believe that this is an ex- 
haustive explanation when both the cruise velocity and 
the elongation of the comets are perpendicular to the 
flat faces of the disks. One must add to the picture the 
links connecting the gel to the host [1]; such links ter- 
minate on the surface at nucleation points , whose dis- 
tribution is determined by the combined action of vari- 
ous activator groups The oversimplified digest given 
above actually translates into a rather intricate dynami- 
cal problem,where the global elastic action of the gel, the 
gel-host links and the polymerization rate come together 



into play. This problem is formulated here in terms of 
a one-dimensional model where the system configuration 
is described by the gel concentration. The gel velocity is 
self-consistently determined by generalizing to the time- 
dependent case the arguments presented in Ref . . The 
time evolution of the gel concentration is described by 
a deterministic equation, which includes elasticity and 
polymerization/depolymerization; the approach being of 
mean-field type, stochasticity is modeled by a diffusion 
term. In our formulation the interaction with the host 
turns into the boundary condition for the concentration: 
the evolution at the boundary has the form of the teth- 
ered ratchet model of Mogilner and Olster jsj. One of 
the main results of the paper is an explicit formula for 
the velocity (see below Eq. [71), wherefrom its correlation 
with the comet size becomes transparent. Furthermore, 
the origin of the threshold mechanism at the onset of 
motion is made clear, as a natural outcome of the model. 
The numerical results shown here exhibit this threshold 
mechanism as well as an intermittent behavior of the ve- 
locity at steady state. The system starts moving past a 
given "latency" time, which depends on the rate of ac- 
cretion of the gel, the elasticity and the nucleation rate. 
The interface polymerization rate can have some rele- 
vance in the initial stages after start, but at steady state 
the bulk contribution prevails: the resulting motion has a 
pulsating behavior, with sudden reductions of the veloc- 
ity and consequent recoveries, as in a system alternating 
"loading" phases at low speed with "unloading" phases 
at higher speed, (see below Fig. [1]). 



II. MODEL 

As anticipated, we are going to establish a one- 
dimensional minimal model for the dynamics of the actin 
gel network during its propulsive action over some host 
object. We take a macroscopic point of view, where the 
concentration and the local velocity are a natural choice 
as dynamical variables describing the gel. Under fairly 
reasonable assumptions one can get rid of the space de- 
pendence of the velocity. The resulting "average" veloc- 
ity is determined by the concentration: hence the system 
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is described by this single variable. In turn the average 
velocity, which we identify with the velocity of the host, 
appears in the evolution equation for the concentration 
(see below EqlH]) . Let us first show how the host velocity 
is determined. When the concentration c{x, t) has a well 
defined peak, say at Xmax, one can approximate the aver- 
age of the velocity field v{x,t) as: < v{t) >« v{xmax,t). 
This obvious remark is not a real progress, unless one 
relates the unknown quantity v{xraax,t) with the value 
of the velocity field at some other point xq where it can 
in principle be determined: such a point is the contact 
of the gel with the host Indeed on the basis of the 
standard rate theories the velocity at the interface is 
a function of the elementary work AM^ according with 
the formula: 



t;A = f+exp(-— -) - (1) 

where V- are the free polimeriza- 

tion/depolymerization velocities. One has 

then: ti(a;o, i) = VA^t)] the link between v{xQ,t) and the 
unknown value v{Xmax,t) is provided by the law of mass 
conservation. From now on, unless otherwise stated, we 
are going to use a reference frame comoving with the 
host, say at velocity v; the mass conservation reads: 

j^cix,t) + -^[{vix,t)~v)cix,t)]=0. (2) 

The role of symmetry breaking in the onset of motion is 
not examined here; let us then assume, e.g., that the gel 
has grown to the left of the host, so that when the host 
moves xq drifts rightwards. We integrate Eql2] over the 
interval Xb < x < xq^ where Xb is a generic position in 
the gel, and obtain: 

r° d 

/ dx—c{x, t) + [v{xo, t)-v]c{xo,t) = [v{xb, t)-v]c{xb, t) 

J Xb 

(3) 

Let us finally assume that the gel and the host drift 
rigidly together, i.e. that v = v{t) = v{xmax,t); upon 
chosing in Eq. [3] the generic position at the peak of the 
concentration (xf, = x^ax) we get the relation we were 
searching for: 

v{t)^v(xQ,t) + - -/ dx—c{x,t) (4) 

This formula gives the velocity v(t) as a nonlocal function 
of the concentration: the evolution equation for c{x, t) 
will make our treatment self-consistent. Let us now dis- 
cuss how such an equation can be determined. The dy- 
namical response of a normal gel, where the structure of 
the network is fixed by the distribution of crosslinks and 
of flexible chains, is characterized by well defined values 
of the coupling constants, such as, e.g., the Young mod- 
ulus in the case of elasticity. A network of semiflexible, 
rather than flexible, chains is apparently more apt to de- 
scribe the actin filaments; whenever its structure is fixed 



its elastic response is also defined, but differs from that 
of a normal gel [IQ] . Apart from this difference one must 
be aware of the fact that the actin network of the comets 
has no predetermined structure: its couplings evolve with 
the configuration, e.g. the elastic response changes as 
the mass and crosslink densities build up or decay: it 
would be appropriate to consider both populations. In 
the present "minimal" context, where we limit ourselves 
to the mass density, the compression modulus depends on 
c(x, t) and the role of crosslinks is included in an effec- 
tive coupling constant C; we refer to the Section Elastic 
Energy for the explicit form of the elastic term. 

The environmental factors, such as the "feeding" of 
polymerization, the depolymerization, the intrinsic fluc- 
tuations of the concentration at the microscopic scale, 
can be assumed to be proportional to the local concen- 
tration. Sure enough, the environment evolves with the 
structure as well, but we will disregard its time evolu- 
tion for the time being. We consider a potential U{x) 
(a growth rate), describing the "feeding" and "decaying" 
processes, in dependence with its sign. Furthermore, a 
diffusive term accounts for the small scale disorder of 
the environment, which perturbs a purely deterministic 
buildup of the network. We add a quadratic term de- 
scribing the short range self-repulsion between different 
portions of the structure. 

Under such factors, in a static reference frame the con- 
centration would then evolve according to the equation: 

d d^ 

c^U(x-vt)c-bc^ + D—c (5) 
dt d-'x 

where the time dependent potential accounts for the mo- 
tion of the " activators" which are are predominantly con- 
centrated around the host. If, e.g., the potential has a 
constant value U > 0. in the feeding region, its action 
combined with the repulsive term —bc^ drives c(x,t) to- 
wards the value c = U /b. The host enforces the condition 
c{x,t) = 0. (cc > Xq); furthermore, the gel is attached to 
it at the nucleation points and grows against it when 
the actin filaments polymerize. We model this along the 
lines proposed by Mogilner and Oster Q, i.e. by means 
of two concentrations a{t),w{t) respectively describing 
the "attached" and" working" portions of the gel. Their 
evolution determines the boundary condition of c at the 
interface:c(a:o, = a{t) + w{t). We now go over and 
write the evolution equation including the elastic term; 
this term is proportional to the deviation of the concen- 
tration with respect to a reference "background" config- 
uration [Tl| . and is opposite to it. If, e .g., an excess 
of mass has accumulated in some region during the pro- 
cess of buildup of the network, the elastic force acts as 
in a spring compressed in the x direction, and elongates 
the excess of mass along x (see the Section Elastic En- 
ergy) . The frame comoving with the host has coordinates 
x' — X — vt, t' = t; in these variables the evolution equa- 
tion has the following form: 
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-U{x)c{x,t) - bc^{x,t) 



6c{x, t) 



(6) 



where we have dropped the primes; Ef.i is the elastic en- 
ergy and r is the relative decay rate. Notice that a drift 
term has appeared, in going to the moving frame: it has 
the effect of elongating the profile of c{x, t) in the direc- 
tion opposite to as if the concentration experienced 
a backward wind. The higher the velocity, the stronger 
the effect: the elongation of the profile "measures" the 
velocity. The velocity and the elongation of the profile 
actually have a link which goes beyond this purely hydro- 
dynamic effect. In fact, according with EqUthe velocity 
V is in turn c-dependent, and upon inserting the l.h.s. of 
EqlB]into Eq[3]we finally get: 



_ v{xo)c{xq) 
^{y^max ) 

dy[U{y)c{y) 



(^\^max ) 



(7) 



where Xmax < xq and the time dependence of the argu- 
ments is understood. In the time independent case only 
the first term survives, and the formula simply states the 
conservation of stationary fluxes [7.] ; the boundary term 
D{j^)xo is negative, as one easily verifies: it averages 
over the region x^ax < x < xq the loss in propulsion 
originating from network's disorder. 



III. NUMERICAL RESULTS 

The formula [7] is one of the primary results of the pa- 
per; it shows how the velocity gets contributions both 
from the interface and from the bulk. Notice that v goes 
as the inverse of the maximal concentration; the elas- 
tic force must counter the formation of exceedingly high 
concentration peaks in order to have higher velocities. 
The bulk contribution becomes increasingly relevant as 
the peak at Xmax moves off xq in the backward direc- 
tion. In the feeding region, close to the interface, we 
have U{x) = U > 0.: the network grows and contributes 
positively to the propulsion. The diffusion constant D 
[D = 5.10~^(/im)^/sec) has a minor effect: its value has 
been chosen two orders of magnitude smaller than the 
diffusion constant of monomers within the gel estimated 
in Ref . . The elastic and self- repulsion terms oper- 
ate positively or negatively dependending on the config- 
uration. They give rise to an intermittent behavior of 
the velocity even at steady state; typically the velocity 
'jumps' between two (or more) values, staying most of 
the time in one of them, which can be identified with the 
cruise velocity of the host. Among the temporary, in- 
termittent velocities we found negative values as well: 
this happens, e.g., when the rate F is half the value 
r = 2. {pN ■ sec)~^ {nm)~'^ used in the simulation pre- 
sented here. At steady state the motion can be depicted 
as a sequence of very short "charging" steps, where the 




FIG. 1: Typical behavior of the velocity (in units of 10. 
nm/sec) versus time (in seconds); notice the discontinuity at 
the onset of motion, the transient before reaching steady state 
and the intermittent behavior at steady state. The values of 
the parameters are reported in the text 



velocity is small or negative, followed by longer advance- 
ment steps. A typical behavior is shown in Fig. [Uwhere 
U = l.(sec)"^: the cruise velocity is f « 8. nm/sec, and 
on average xq — Xmax — 250. nm. The velocity is quite 
near to 10. nm/sec, established in Ref.[l| as a minimum 
value for classifying an object as being in motion. At 
a closer inspection the concentration profile reveals time 
dependent small size deformations tuned with the steps of 
the velocity. In particular the peak position Xmax evolves 
with an intermittent behavior: temporary compressions 
are subsequently released. Let us call active region the 
portion of space going from the peak to the interface; 
in the example shown in Fig[T] this region goes slightly 
beyond the feeding region, which has size 200. nm (see 
FigE]). The formula for the velocity Eq. [7] shows that 
the bulk contribution is proportional to the size of the 
active region. A linear relation between the length of the 
comet and the cruise velocity has been observed in very 
accurate kinetics experiments by Soo and Theriot [l| : we 
believe that our results arc consistent with such observa- 
tions. Actually there are two concomitant effects: the 
drift term elongates the profile by a factor proportional 
to V, in turn a more elongated profile generates a larger 
V. Let us examine the interface contribution. According 
to Eq. [H v(xq) is a decaying function of the elementary 
work AW: 



AVF 



51 ■ F{xo) = SI ■ (Tx.x{xo) 



SI ■ c{xo) — Ca 
w 



dy[c[y) - c{xo)] (8) 



where SI is the elementary displacement having the or- 
der of the size of a monomer SI = 2.2 nm, — the ratio of 
attached versus working gel and Ca = 2. [pN ■ nm) the 
elastic coupling of the attached portion of the gel. The 
above expression follows from the equilibrium condition 
between working forces and load Q, the load here be- 
ing uniquely given by the attached links (see the Section 
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FIG. 2: Latency time (in sees) versus the growth rate U (in 
units (sec)~^);the size of the feeding region is Xf — 200. nm; 
notice that when U <U ~ 0.6 the motion cannot start; these 
results where obtained with F — 0.1 {nm)~^ ■ {pN sec)~^ , Ca ~ 
l.pN ■ nm, Cw = 0.5 pN ■ nm; all the other parameters have 
values as reported in the text 



Elastic Energy). Let us examine the order of magnitude 
of the force at steady state in the case reported in Fig. 
[TJ There the ratio a/w is close to one (see Fig. m, the 
active region has size xq — Xmax — 250. nm, the peak is 
c{xmax) ~ 10. (the concentration has units (10. nm)^^) 
the deviation c{xmax) ~c{xq) is of order 10^^ so that one 
has F{xi^) ~ 5.pN. Taking into account that the Boltz- 
mann factor is fcsT ~ A.lpN ■ nm one can verify that 
the contribution to v arising from v[xo) is quite small as 
compared with the bulk contribution. 

Assuming a free polymerization velocity t)+ = 
10. nm/ sec and w_ — 0. we always obtained that the bulk 
velocity has the dominant role in propulsion, with the ex- 
ception of the initial transient following the onset of the 
motion. Our simulations start with zero concentration. 
A given nucleation rate drives the attached gel (a(t)), its 
decay activates w{t) Q. As a{t) and w{t) start increas- 
ing, the boundary condition c{xo,t) = a{t) + w(t) initi- 
ates the evolution of c{x, t). The velocity does not follow 
smoothly the gel growth, but rather jumps abruptly from 
zero to a finite value: the onset of motion has a threshold. 
The motion can start provided that the gel generates a 
positive propulsion: this happens when the maximal con- 
centration rises above c{xo, t): at this point an excess of 
mass has been created on the back of the host, and the 
elastic response to it translates into motion. If the gel 
growth rate is too small, the maximal concentration stays 
locked at the interface, and the motion cannot start. 

If the growth rate is sufficiently high the maximum of 
c stays at xq, but only for a finite "latency" time. This 
time is strongly dependent on the rate U in the feeding 
region; in Fig.[2]it is shown that it diverges as U decreases 
towards a finite value U « 0.6. Similarly, one can ver- 
ify that shorter feeding regions generate longer latency 
times. As the peak unlocks from xq, the host is instantly 
put into motion, and its velocity v abruptly departs from 



20 60 100 

FIG. 3: From bottom to top: potential U(x); profile c(x) at 
t=2.sec, right before the onset of motion (dotted line); c(x) 
at steady state (t=15. sec) (continuous line). The abscissa has 
space unit a = 10. nm. The vertical axis is in units {sec)~^ 
(U(x)), and (10. nm)-i (c(x)) 




FIG. 4: Time behavior of the variables a (continuous line) and 
w (dotted line), representing the attached and working part 
of the gel at the interface; same parameters as in Fig. [T]the 
vertical axis has units (10. nm)~^, the time unit is l.sec. At 
t=l. w is still sensibly smaller than a. Notice that at some 
time t ^ 2. w grows above a: this event is close to the onset 
of motion 

zero, with a value depending on the mass accumulated 
during the latency. In the evolution equation Eq. [H] the 
drift term is simultaneously activated: the profile is faced 
with an abrupt backward wind, and must adapt its shape 
to it. A transient follows, during which the active region 
modifies its size with strong fluctuations, until the system 
reaches a steady state profile. In Fig. [3] we compare the 
profile during latency, but close to start, with the steady 
state profile; the potential is superimposed \U{x) = 1. in 
the feeding region, U{x) = —0.2 in the decay region]. The 
steady state motion results from the combined action of 
the drift term and of the elastic forces; the drift term 
in particular marks the shape of the steady state profile, 
which appears smoothly elongated away from the inter- 
face as compared with its form at the onset of motion. 
A stationary transfer of mass from the feeding region to 
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the queue maintains the steady state profile. 



IV. ELASTIC ENERGY 

We write the elastic energy as a quadratic form 

Eel = ^ J dx J dy Sc{x, t) ■ G{x, y) ■ Sc{y, t) (9) 

where 5c{x,t) — c{x,t) — c{x,t) is the deviation from a 
"background" configuration c{x, t), C is an energy whose 
value must account of the effect of crosslinks. The ker- 
nel G{x,y) has the role of bounding the range of elastic 
response within a finite length A: G{x, y) — exp( ). 
We define background or reference configuration a config- 
uration which is in static equilibrium but only over a suit- 
able time scale. New configurations generated at shorter 
time scales are attracted towards it. In the present pa- 
per we give a simplified treatment where the elastic in- 
teraction is limited to the active region: we believe that 
this region has the dominant role in the propulsion. In 
the bulk of the network the elasticity operates against 
changes in the concentration profile. The background 
configuration results from the past hystory of the struc- 
ture. For the time being we disregard memory effects 
and simply mimic the rigidity of the system with the 
self-repulsion term included in Equation [BJ the relative 
coupling constant has value b = 0.5 nm/ sec. Let us fo- 
cus now on the active region: there a natural choice for 
the background configuration is given by the concentra- 
tion at the interface. This is because at the host's sur- 
face the concentration is fixed by the boundary condition 
c{xQ,t) = a{t) + w{t), so that in the adjacent region the 
profile must be predominantly attracted towards c{xo, t). 
In specializing the energy Eq. [9]to the active region we let 
then c{x,t) = c(xo,t);it is understood that the integra- 
tions are extended over the interval Xmax < x < xq, and 
we assume A >> (xQ—Xmax)- About the form of the cou- 
pling C, one must consider that the attached and work- 
ing portions operate in opposite ways and with different 
strength. The variable a{t) is related with links connect- 
ing the gel to the host: their role is of the same nature as 
that of the gel-gel crosslinks, which determine the bulk's 
elastic response. Hence a{t) attracts the system towards 
the background configuration while w(t) acts in the oppo- 
site direction. The specific form of the coupling C must 
then be: C = ■^^Ca--^^C^u where with Ca and C^, we 
have separated the elastic response of the attached links 
from that of the working filaments (we always assume 
that Ca > Cw'i in particular in the simulations presented 
here we have Ca = Z.Cu,; Cw = l.pN ■ nm)). As we as- 
sumed that the length A is much larger than the active 
region, we average over the space dependence of the devi- 
ation: Sc{x,t) Ki Sc = ^{c{xmax,t) — c{xo,t)). The total 
elastic force, proportional to Sc, stabilizes the gel pro- 
vided that the contribution from the gel-host links is the 
dominant one, i.e. under the condition Ca -a — Cw - w > Q. 



Let us add some comment on the dynamics of a{t) and 
w{t): their evolution equations are slightly modified with 
respect to equations (1) and (2) of Ref. Here w{t) 
gets a contribution not only from the detached links, but 
also from the adjacent gel (c = c{x)\ {xq — x) « 1): 



i - 5{v)a 
S{v)a — dcW 



D{c~2w) 



(10) 



here 5{v) = 5o{\. + |) is the velocity-dependent detach- 
ment rate, dc is the capping rate and j is the nucle- 
ation rate per space unit. We assigned to these variables 
the following values:j = 0.5 {nm ■ sec)^^ 5q = 1. (sec)~^ 
V = 100. nm/sec, dc = 0.8 (sec)^^. It is necessary to 
connect this formalism with the usual treatment of elas- 
ticity: this is done by assuming that the system is in- 
compressible, i.e. that the elementary mass contained 
in a generic interval Ax is conserved under deformation: 
d{c{x)Ax) = 0. The condition above gives: 



Sc{x) 



-c{x) — Ax 
dx 



-c{x)Ux^xix) 



(11) 



where Ux.xix) is the strain tensor in one dimension. Eq. 
|9]then reads: 



C 

Eel ^ — I dx 



dy[c{x)ua:^x{x)]G{x,y)[c{y)uy^y{y)] 



(12) 

this should be compared with the standard form of the 
elastic energy: 



Eel = B dx{ux.x{x)Y 



(13) 



The compression modulus which in one dimension is 
an energy/length, in our case is explicitly concentration- 
dependent. In fact, the stress tensor ax,x{x)^ which here 
has the dimension of a force, is given by 

5Eel 



^X,X (^) 



5ux,x{x) 

c(a;)-^^2L 
5c{x) 



CXeffC^{x)Ux,x{x) (14) 



where Ae// is the effective range of the interaction. We 
have then: B CXef fc^{x). In the case of Fig [3] one 

would have B W Ca ■ {xq - Xmax) ■ C^iXmax) ~ 0.5 71 A^. 

Correspondingly the elastic stiffness can be estimated as 
the ratio fcgj = B/{xo — Xmax)' we obtain then kei = 
10^^ nN/nm, quite small as compared with the value 
kei — 0.17 nN/nm, reported in Ref. |13j |. 

The value of Ca was chosen by requiring a AW^ of order 
ksT. In the equation of motion the elastic coupling C 
multiplies F, whose value was chosen in order to make 
the elastic term strong enough to counter an unlimited 
growth in the feeding region. 

From EqlS] one can argue that a steady state with a 
ratio a/w << 1. would allow for higher values of Ca and 
consequently would give more realistic values of stiffness; 
the rate F should be accordingly reduced, in order to keep 
the same strength of the elastic term in the equation of 
motion. 
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V. CONCLUSIONS 

We presented a continuum model of actin-mediated 
bacterial propagation. The model is a minimal one, in 
that it describes the system in terms of a single variable: 
the concentration c{x,t). The time evolution of c{x,t) 
includes the interaction with the environment and with 
the host as well as the elastic interaction in the region 
adjacent to the host's surface. A formula for the velocity 
is obtained, where the interface and bulk contribu- 
tions are made explicit. This formula, together with the 
hydrodynamic effect acting on the profile, makes the cor- 
relation between the velocity and the length of the actin 
network explicit: this result agrees with the behavior 
observed in Ref. and indicates that the bulk, rather 
than the interface, has the major role in propulsion. The 
model furthermore motivates on theoretical grounds the 
existence of an intrinsic threshold mechanism underlying 



the onset of the motion. A relevant feature of a minimal 
model is that in the numerical simulations the number 
of adjustable parameters is accordingly quite small. 
The key parameters regulating the dynamics of the 
concentration in the bulk are the growth rate f/,the 
self-repulsion coupling b, the elastic coupling C and 
the rate F; the interface dynamics is controlled by the 
nucleation, detachment and capping rates (j, S{v), dc). 
Preliminary numerical results exhibit the threshold 
mechanism and its dependence on the growth rate (see 
FigsUl [2|. The "cruising" steady state is characterized 
by a pulsating behavior, where the motion undergoes 
sudden reductions and subsequent recoveries. 
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